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Abstract 

Using a simple model with link removals as well as link additions, we show that an evolving network 
is scale free with a degree exponent in the range of (2,4]. We then establish a relation between the 
network evolution and a set of non-homogeneous birth-and-death processes, and, with which, we capture 
the process by which the network connectivity evolves. We develop an effective algorithm to compute the 
network degree distribution accurately. Comparing analytical and numerical results with simulation, we 
identify some interesting network properties and verify the effectiveness of our method. 
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Growing Networks and Pure Birth Processes 

The first growing network model, i.e., the BA model proposed by Albert, Barabadsi, Jeong 
Q], predicts a power-law network degree distribution with exponent 7 = 3, whereas the degree 
exponents of many real complex networks are found empirically in the range of (2, 4) . This 
has motivated extensive research to modify the basic BA model to match with practical scale-free 
networks. These variations can be summarized by the following general BA model: 

(i) Initialization: n Q nodes are given at time t — 0; 

(ii) Growth: at the tth time step, a new node and m(t)(< t + no) new links from this node are 
added; 

(iii) Preferential attachment: the new node is connected to an existing node i according to the 
following probability H(ki) = f(ki)/ f(kj), where ki is the number of degrees of node i and 

is pre-selected function. 

With general m(t) and f(ki), this model is analytically intractable and is very complicated for 
simulation study. We may simplify this model in two ways. Setting m(t) = m, i.e., the number of 
links in the network increase linearly, the growth of the network is stationary. Assuming f(ki) = ki 
further, we have the basic BA model. If we set f(ki) = (1 — p)ki + p instead, where p is the 
probability of randomly selecting an old node, the model reduces to Liu et al. jjj's model with 
7 = 3 + p/[m(l — p)]. In Bianconi and Barabsi p's fitness model, f(ki) = rjiki, where r\i is 
chosen from a distribution p(rj). If p{rf) is uniform, 7 = 2.255. Alternatively, we can also first set 
f(ki) = ki. With a time-dependent number of new links added at each time step, the growth of 
a network is non-stationary. For example, with the accelerating function m(t) = mt e , < 9 < 1 
proposed by Dorogovtseva and Mendes (|, the degree exponent 7 = (3 — d)/(l — 9) and the non- 
stationary exponent z = 29/(1 — 9). Shi, Chen and Liu \L\ proposed a slower accelerating function 
m(t) = mint, t > 2. 

n 

Shi, Chen and Liu |7( established a relation between the connectively of a growing network and 
a set of non-homogeneous pure birth processes (PBP) and found numerically that 7 ~ 3.1 and the 
non-stationary exponent is very small for the case of m(t) = mint. It can be observed that the 
deg.ee distribution curves of growing networks inB is snake Irke with a slightly downward headrag 
head section as illustrated in figure [T] (see, also |2| for similar figures). 

There are other ways to extend the BA model. Readers can refer to jsj for a comprehensive 
review. 
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FIG. 1: The degree distribution for the BA model with m — 5, S — 3 and t = 150000: by the analytical method in 
solid green line; by simulation in dotted black line; and by the PBP method in solid red line. 



A Model of Evolving Networks and Dynamic Equation 

Albert and Barabadsi js| considered a model of the evolving network in which some old links are 
rewired at each time step (see Section 4). We observe from many real networks that beside adding 
new nodes and links, some old nodes and links can also be removed as a network evolves. In other 
words, many networks display a dynamic evolving process. We propose the following simple model 
to capture the basic features of the above evolving network. 

(i) Initialization: There are n fully connected initial nodes. 

(ii) Link removal: At each time step, c old links are removed as follows. We first select node i 



with the anti-preferential probability similar to that used in |2| 



The link connecting nodes i and j is removed. We repeat this procedure c times to remove c existing 
links. Finally, isolated nodes are removed from the network. 

(iii) Link addition: At each time step, a new node is added to the system and m(< no) new links 
from the new node are connected to m different existing nodes. A node i with degree ki will receive 
a connection from the new node with a Bayes' preferential probability 





where a is used normalized factor such that a 



-i 



J^. k i . We then choose node j from the 



neighborhood of node i (denoted by Oi ) with probability K i 1 Y\*{kj) 1 where Ki = YljeOi 



n(fci) 



fci + i 



(2) 
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The above model is different from the existing ones in that it removes links instead of rewiring 
links. Furthermore, isolated nodes are removed in our model. 

By the continuum theory, ki(t) approximately satisfies the following dynamic equation 



— = mU( kl 



h + l 2 

~ mi—, ; ; c-, (6) 

[2(m - c) + l]i t W 

where the last approximation is based on ^2j{kj + 1) = 2(m — c)t + N(t — 1) ~ [2(m — c) + l]t, 
J2jt=Oi -^7 1 n*(^j) ~ 1, and, in the mean-field sense, ak^ 1 « [N(t - « 1/t, in which N(t - 1) 
is the number of non-isolated nodes at time step t. 

Let ti be the time step when node i is added to the network. Initially, node i has /^(tj) = m 
links, thus the above equation has the following solution 

Xi (t) = \ki(t) + B-m\ =B (£\ , (4) 



with the dynamic exponent 



and the dynamic coefficient 



TYl 



m - 2c 2(m - c) + 1 
B = B(m,c) =m + 1 -± >— -i 6 

m 

In the solution procedure, we require < /? < 1 and _B > for the solution to be feasible. Some 
simple analysis of the above formulas shows that m > 2c is a sufficient condition for equation © 
to have a feasible solution. 

Assume that ti follows a uniform distribution over interval (0,t). We have, by 

P(x) = ^-B 1/p x-\ (7) 

P 

where x G [B, oo) following from (pEJ) and the degree exponent 

7 = 1 + ^ = 3 + ^- (8) 

p m 

Equation (jSJ) shows that this system self-organizes into a scale-free network with 2 < 7 < 4. 
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The next step is to obtain the network degree distribution. For B > m, we have, following the 
standard mean field approach |8[, the explicit solution of dynamic equation (jSJ) 



ki{t) = B 




1 



+ m 



(9) 



and the network degree distribution 



P(k) = -B^^k + B -m)~ 7 . 



(10) 



For B < m, the continuum theory does not render an accurate solution, and we need a different 
method. 

Birth-and-Death Processes of Network Connectively 

The dynamics of the degree of a node in an evolving network is closely related to Markov 
processes. Let K,i(t) be the degrees of node i at time t. Since Ki(t + 1) only depends on Ki(t) 
and allows the removal of old links for our model, {Ki(t),t = i,i + 1, . . .} is a discrete-time Markov 
process with the state space Q = {0, 1, 2, . . .}. 

By (jnj), the probability that node i with degree k(> 1) is connected to a new node at time step t is 
approximately gt{k) = (l~2c/t)mll(k) « mll(k) = m(k+l)/([2(m — c) + l]t). The probability that 
note z's degree decreases by 1 is approximately (2c/ 1) [1 — mU(ki)} rs 2c/ 1, while the probability 
that its degree decreases by more than 1 is o(t) and will be ignored. Thus, the probability that 
the degree of node % remains the same is h t (k) = 1 — g t (k) — 2c/ 1. This shows that K^(t) is in fact 
a non-homogeneous birth-and-death process (BDP). In addition, we set poo = 1 since we remove 
isolated nodes and pkk = 1 when k > m + t — i. In summary, for t — i, i + 1, i + 2, . . ., the one-step 
transition probability matrix of node i at time t is given by 



1 







2c/t ht(l) 



9 t (l) 



Pi(* + 1) 



(11) 



2c/t h t (m + t — i) g t (m + t — i) 







1 
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Denote = P{Ki(t) = n} for n = 0,1,2,... and = (/ i(0 (t), /;,i(t), • • • , / 4 , n (t), • • •)■ 

Obviously, /^(i) = (0, 0, • ■ ■ , 0, 1, 0, • • ■) = e m where ej jm (i) = 1. By density evolution, the (t + l)th- 
step probability vector fi(t + 1) for node i is given by 

,A(t+l)=e^-P i (i + l)-P < (i + 2)...P i (t+l), t = i,i+l,.... (12) 

Let 

p{s,t)( t + i) = J2 fi(t + 1) = g rrPi{i + 1) • P«(i + 2) • • • P<(t + 1), (13) 

where the integer 5 > 1 is needed technically for the transition probability matrix. For the choice 
of S and its impact on computation, please refer to 

Generally, it would be extremely difficult to calculate (|T3j) . Fortunately, we can find the following 
relations: 

e m Pi(t) = e m P s (t), i = S + l,S + 2,...- t = i + l,i + 2, (14) 
and, in general, for s = 1, 2, ... 

e m Pi(t)Pi(t + 1) • • • Pi(t + s) = e m P s (t)P s (t + 1) • • • Px(t + s). (15) 

Thus we obtain the following key algorithm 

F(W(t + 1) = ((• ■ ie m P s (S + 1) + e m )P s (S + 2) + ■■■) + e m ) ■■ ■ P s (t + 1). (16) 

The right-hand side of (|16j) can be efficiently computed with a complexity of 0(t 2 ) [7(. 

The degree distribution of a network can be determined by the average of the degree distributions 
of all the nodes. Therefore, for a sufficiently large t, we have 

P(k)*P(k,t + l)= f^fZf- . (17) 

As a bonus, we can also estimate the number of non-isolated nodes from (fTTj) as follows 

N(t)&(t + l)[l-P(0,t + l)}, (18) 

noting that P(0, t + 1) is the probability that a node is isolated at time step t+1. This index cannot 
be obtained from (fTtjjh 

We note that there are also inaccuracies in the transition probability matrices, and we can only 
perform a finite number of computation steps to estimate the asymptotic network behavior. To 
verify the BDP method, we compare the computation results with simulation. 
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FIG. 2: The degree distribution obtained by (|10|) in solid green line; by (|19|) in dash blue line; by simulation in 
dotted black line; and by the BDP method in solid red line for m = 5, c — 1, S — A and t = 150000. The settings for 
the inserts are the same as those of the main figure, except m = 3 and c = 1 for the bottom left insert and m = 8 
and c — 2 for the top right insert, and no simulation results are given. 

From figure 121 we see that the network degree distribution curves obtained by the BDP method 
and by simulation match very well. Our method predicts a horse head-like distribution curve, with 
its middle section displaying the expected scale-free state. Because we can only perform a finite 
number of computation steps, there is an inward bend at the tail of the distribution curve (see j]| 
for more detailed discussion). The degree exponent and coefficient can be estimated by applying 
the least square method to the data generated from P(k) G (10 -4 , 10~ 6 ). 

The numerical results from the birth-death processes clearly show that the network degree dis- 
tribution curve has a very different head section. This motivates us to construct the following 
approximation for the degree distribution when B < m, noticing also from (@J that, when B < m, 
Xi(t) achieves the minimum at = m and is symmetric 




k > m 



k < m 



(19) 



where /i is a fitted parameter and the coefficient 




(20) 



is a normalizing constant such that J °° P(k)dk = 1. 

Empirically, we find that when u = 0.2m + (c — 1), the approximation (fT^j) is very accurate for 
the overall distribution and captures the pattern of the small degree distribution, as shown in the 
small inserts in figure 2. The figures also show that (|10|) cannot provide probabilities for degrees 
smaller than m, and it visibly over estimates other small degree probabilities and under estimate 
large degree probabilities. 

Application to the Albert-Barabdsi model 

The model proposed by Albert Barabdsi |a] starts with no isolated nodes, and performs one of 
the following operations at each time step: 

(i) Add m(< no) new links with probability p: Select a node randomly as the starting point of 
the new link and then select the other end of the link with the preferential probability (0). Repeat 
this process m times. 

(ii) Rewire m links with probability q: Select randomly a node i and a link connected to it. 
Remove this link and replace it with a new link l^i that connects % to node j' which is chosen with 
the preferential probability (0). Repeat this process m times. 

(iii) Add one new node with probability r = 1 — p — q: The new node has m new links that are 
connected to different existing nodes with the preferential probability (J2J). 

By the continuum theory, Albert and Barabdsi obtained the following dynamic equation 

dh ki + 1 p-q , 

— — ~ m- ; r : — V m , (21) 

dt [2m(l -q)+r]t rt K J 

and from which they derived the network degree distribution 

P(k) = -{m + T)?{k + T)-\ (22) 

where r = (p — q)(2m(l — q)/r + 1) + 1 and the degree exponent 

1 r — 2mq , , 

7 = 1 + -, = 3 + q -. 23 

p m 

Obviously, (}2*2"j) is valid only when m + r > 0. Thus for the Albert-Barabadsi model, the network 
degree distribution is scale-free only when parameters p and q satisfy 

q < q m ax = min{l - p, (m + 1 - p)j (2m + 1)}. (24) 
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Now, we use birth-and-death processes to discuss the Albert-Barabadsi model. By (|21|) . the 
one-step transition probability matrix of node i at time t is given by 

ri-<fc(o) &(o) i 



mq/rt h t (l) 



9t(l) 



Pi(* + 1) 



(25) 



mq/rt ht(m + t — i) g t (m + t — i) 







1 







where gt{k) ~ (m(k + 1) /([2m(l — q) + r]t)) + mp/rt and /i t (fc) = 1 — gt(/c) — mq/rt. 

The results from (J2~2"|) and the BDP method are compared in figures |3] and 0J The distribution 
curves of the two methods match very well in the middle section. Again, the distribution curves 
from the BDP method are horse head like with a downward bending head. (j22J) over estimates small 
degree probabilities and does not provide the probabilities for degrees smaller than m. 



FIG. 3: The degree distribution obtained by (|22|l in solid green line and by the BDP method in solid red line for 
m = 2,p = 0.6, q = 0.1, S = 2 and t = 100000. 
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FIG. 4: The degree distribution obtained by (|22|l in solid green line and by the BDP method in solid red line for 
to = 5, p = 0.2, q = 0.4, S = 4 and t = 150000. 

We summarize the results and findings in this paper as follows: (1) We introduce a simple yet 
flexible model of evolving networks with both addition and removal of links and nodes. The removal 
of both links and isolated nodes is new; (2) The connection between an evolving network and a 
set of non-homogenous birth-and-death processes provides an efficient algorithm to numerically 
calculate the network degree distribution. With this method, we reveal the complete process by 
which a network evolves into a scale-free state; (3) With the close match between the numerical 
results and simulation results, our birth-death method provides an efficient and reliable substitution 
to simulation, in particular since the existing analytical methods cannot handle more complicated 
network mechanisms and the computational requirements of simulation are often too high; (4) We 
find that the method based on the continuum theory is not suitable for small degree distribution 
and under estimates large degree probabilities; (5) The horse head-like degree distribution curves 
have been observed in a number of real networks, such as the actor collaborations and word co- 
occurrences networks (see figure 1 (d) and (e) in Newman Using the birth-and-death process 
method, we demonstrate that the distribution curves of growing network are snake head like while 
the distribution curves of evolving networks are horse head like; and (6) Degree distributions of 
evolving networks have two distinct sections and the maximum probability occurs at degree m. 
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